System and method for single-beam internal reflection tomography

ABSTRACT

A methodology and concomitant system for three-dimensional near-field microscopy achieves subwavelength resolution of an object via total internal reflection. The features of this approach include: (i) the evanescent waves used for illumination encode on the scattered field the subwavelength structure of the scattering object—it is thus possible to obtain subwavelength resolved images of the sample as is done in other near-field techniques such as near-field scanning optical microscopy without the technical difficulties encountered with probe-sample interactions; and (ii) the results of the reconstruction are unambiguous in the sense that the relation between the scattered field and the three-dimensional structure of the sample, as described by the spatial dependence of the susceptibility, is made manifest.

BACKGROUND OF THE DISCLOSURE

1. Field of the Invention

This invention relates to tomography and, more particularly, totomography wherein an image of an object is directly reconstructed withsub-wavelength resolution.

2. Description of the Background Art

Total internal reflection microscopy (TIRM) brings to conventionalmicroscopy the added functionality of illumination by evanescent waves.Incorporation of evanescent waves in the illuminating field is animportant development for several reasons. First, the exponential decayof such waves along one direction allows for the control of the depth ofpenetration of the illuminating field. Second, evanescent waves may beused to resonantly excite surface plasmon modes of the sample. Finally,and perhaps most importantly, evanescent waves may be employed tosupercede the Rayleigh diffraction limit of order half the wavelength,λ/2.

TIRM has been in practical use for decades. It has primarily been usedas a surface inspection technique, though the sensitivity of the fieldto distance along the decay axis has been utilized to advantage inapplications such as the measurement of distance between two surfaces.Until recently the opportunities for transverse superresolution madepossible by the high spatial frequency content of the probe field havebeen largely over-looked. However, there have been advances in therealization of super-resolved imaging utilizing TIRM. There had been aninitial direct imaging approach resulting from the marriage ofstanding-wave illumination techniques and TIRM achieving transverseresolution of λ/7. A later approach in which the information content ofthe TIRM experiment for the case of weakly-scattered fields has beendetermined is reported by Fischer in “The information content of weaklyscattered fields: implications for near-field imaging ofthree-dimensional structures”, published in Journal of Modern Optics,Vol. 47, No. 8, pages 1359-1374, 2000.

However, the art does not teach or suggest a technique for thereconstruction of a tomographic image of the object which isaccomplished by making use of an analytic solution to an inversescattering formulation with evanescent waves.

SUMMARY OF THE INVENTION

These shortcomings, as well as other limitations and deficiencies, areobviated in accordance with the present invention via a technique,referred to as total internal reflection tomography (TIRT,) whereby anobject is probed with evanesent waves and then waves scattered from theobject are processed with a prescribed mathematical algorithm toreconstruct the tomographic image.

In accordance with a broad method aspect of the present invention, amethod for generating a tomographic image of an object includes: (1)probing the object with incident evanescent waves; (2) detecting wavesscattered by the object, and (3) reconstructing the tomographic image ofthe object by executing a prescribed mathematical algorithm withreference to the incident evanescent waves and the scattered waves togenerate the tomographic image with subwavelength spatial resolution.

A broad system aspect of the present invention is commensurate with thebroad method aspect.

The features of the TIRT approach are at least three-fold: (i) theevanescent waves used for illumination encode on the scattered field thesubwavelength structure of the scattering object. It is thus possible toobtain subwavelength resolved images of the sample as is done in othernear-field techniques such as near-field scanning optical microscopywithout the technical difficulties encountered with probe-sampleinteractions; and (ii) the results of the reconstruction are unambiguousin the sense that the relation between the scattered field and thethree-dimensional structure of the sample, as described by the spatialdependence of the susceptibility, is made manifest. This is somewhatanalogous to the transition from projection radiography to computedtomography.

BRIEF DESCRIPTION OF THE DRAWING

FIG. 1 illustrates a object/scatterer under test showing inputevanescent waves generated at a prism face by total internal reflectionan whereby the object scatters evanescent modes to the far zone;

FIGS. 2A and 2B depict direct reconstruction results for an exemplaryscatterer using different parameters;

FIG. 3 is a high-level block diagram of a system for directlyreconstructing the tomographic image of the scatterer; and

FIG. 4 is a flow diagram of the methodology for directly reconstructingthe tomographic image of the scatterer.

DETAILED DESCRIPTION 1. Forward Problem

In this section we consider the scattering of evanescent waves fromweakly-scattering dielectric media. In order to develop the physicalideas in a simple setting, we begin with the case of scalar waves. Wethen turn to the full vector theory of electromagnetic scatteringthereby accounting for the effects of polarization. Note that the vectortheory is essential since the scalar approximation to the scattering ofelectromagnetic waves is invalid when the dielectric susceptibilityvaries on subwavelength scales.

1.1 Scalar Case

Begin by considering an experiment in which a monochromatic scalar fieldis incident on a dielectric medium with susceptibility 72(r). The fieldincident on the sample will be taken to be an evanescent wave which isgenerated by total internal reflection at the interface of twohalf-spaces. One half-space, taken to be z≧0, will have the vacuum indexof refraction 1 while the z<0 half-space will have an index ofrefraction n. The scalar field U(r) obeys the reduced wave equation

∇² U(r)+k ₀ ² n ²(z)U(r)=−4πk ₀ ²η(r)U(r),  (1)

where k₀ is the free-space wavenumber, n(z) is the z-dependent index ofrefraction as described above, and the support of η(r) is contained inthe z≧0 half-space. We represent the total field as the sum

U=U _(i) +U _(s),  (2)

where U_(i) and U_(s) represent the incident and scattered fields,respectively. The incident field obeys the homogeneous equation

∇² U _(i)(r)+k ₀ ² n ²(z)U _(i)(r)=0.  (3)

The scattered field obeys the equation

∇² U _(s)(r)+k ₀ ² n ²(z)U _(s)(r)=4πk ₀ ²η(r)U(r).  (4)

Equation (4) may be recast as the integral equation

U _(s)(r)=k ₀ ² ∫d ³ r′G(r,r′)U(r′)η(r′),  (5)

where G(r,r′) is the Green's function.

The Green's function G(r,r′) satisfies the equation

∇² G(r,r′)+n ²(z)k ₀ ² G(r,r′)=−4πδ(r−r′)  (6)

and obeys the boundary conditions

G(r,r′)|_(z=0) _(⁺) =G(r,r′)|_(z=0) _(⁻) ,  (7)

{circumflex over (z)}·∇G(r,r′)|_(z=0) _(⁺⁼) {circumflex over(z)}·∇G(r,r′)|_(z=0) _(⁻) .  (8)

It may be seen that G(r,r′) admits the plane-wave decomposition$\begin{matrix}{{{G( {r,r^{\prime}} )} = {\frac{i}{2\quad \pi}{\int\quad {{^{2}{{qk}_{z}^{- 1}(q)}}\{ {1 + {{R_{1}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \} {\exp \lbrack {{{ik}(q)} \cdot ( {r - r^{\prime}} )} \rbrack}}}}},} & (9)\end{matrix}$

where R₁(q) is the reflection coefficient given by $\begin{matrix}{{R_{1} = \frac{{k_{z}(q)} - {k_{z}^{\prime}(q)}}{{k_{z}(q)} + {k_{z}^{\prime}(q)}}},{with}} & (10) \\{{k_{z}(q)} = \sqrt{k_{0}^{2} - q^{2}}} & (11) \\{{{k_{z}^{\prime}(q)} = \sqrt{{n^{2}k_{0}^{2}} - q^{2}}},} & (12)\end{matrix}$

and k(q)=(q, k_(z)(q)). The plane wave modes appearing in Equation (7)are labeled by the transverse part of the wavevector q. The modes forwhich |q|≦k₀ correspond to propagating waves while the modes with |q|>k₀correspond to evanescent waves.

Consider now the situation in which the incident field is an evanescentplanewave of the form

U _(i)(r)=exp[ik(q)·r],  (13)

where q labels the transverse part of the incident wavevector. Withreference to FIG. 1, if incoming wave 101 traverses prism 105 with indexof refraction n, then k₀≦|q|≦nk₀. Thus k_(z) is imaginary with thechoice of sign dictated by the physical requirement that the field decayexponentially with increasing values of z. In the first Bornapproximation we replace the field by the incident field in the righthand side of equation (4) and obtain the following expression for thescattered field, which occurs because object 100 frustrates totalinternal reflection and scatters evanescent modes to the homogenous modethat propagates to the far field:

U _(s)(r)=k ₀ ² ∫d ³ r′G(r,r′)U _(i)(r)η(r′).  (14)

In the far zone of the scatterer, the leading term in the asymptoticexpansion of the Green's function is given by the expression$\begin{matrix}{{{G( {r,r^{\prime}} )} \sim {\{ {1 + {{R_{1}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \} \frac{\exp ( {{ik}_{0}r} )}{r}{\exp \lbrack {{- {{ik}(q)}} \cdot r^{\prime}} \rbrack}}},} & (15)\end{matrix}$

where k(q) is parallel to r and |q|≦k₀. Using this result and equation(14) we find that the scattered field behaves as an outgoing homogeneouswave of the form $\begin{matrix}{{{U_{s}(r)} \sim {\frac{\exp ( {{ik}_{0}r} )}{r}{A( {q_{1},q_{2}} )}}},} & (16)\end{matrix}$

where q₁ and q₂ are the transverse parts of the incident and outgoingwavevectors, respectively. Here A(q₁,q₂), which is the scatteringamplitude associated with the scattering of evanescent plane waves withtransverse wavevector q₁ into homogeneous plane waves with transversewavevector q₂, is related to the susceptibility of the scatting objectby the expression

A(q ₁ ,q ₂)=k ₀ ² ∫d ³ r{1+R ₁(q ₂)exp[2ik _(z)(q ₂)z]}exp{i[k(q ₁)−k(q₂)]·r}η(r).  (17)

We will take up the inversion of this integral equation in Section 2.1.

1.2 Vector Case

We now turn to the vector theory of electromagnetic scattering. Asbefore, we consider an evanescent wave which is incident on a dielectricmedium with susceptibility η(r) and assume that the evanescent wave isproduced at the interface of two half-spaces. We consider onlynonmagnetic materials and accordingly will limit our attention to theelectric field E which satisfies the reduced wave equation

∇×∇×E(r)−k ₀ ² n ²(z)E(r)=4πk ₀ ²η(r)E(r),  (18)

where k₀is the free space wave number and n(z) is the background indexof refraction as described in Section 1.1. We decompose the field into asum of two parts,

E=E ^(i) +E ^(s).  (19)

The incident field E^(i) obeys the homogeneous equation

∇×∇×E ^(i)(r)−k ₀ ² n ²(z)E ^(i)(r)=0.  (20)

The scattered field E^(s) obeys

∇×∇×E ^(s)(r)−k ₀ ² n ²(z)E ^(s)(r)=4πk ₀ ²η(r)E(r).  (21)

Equation (21) may be reformulated as the integral equation

E _(α) ^(s)(r)=k ₀ ² ∫d ³ r′G _(αβ)(r,r′)E _(β () r′)η(r),  (22)

where G_(αβ)(r,r′) is the Green's tensor for the half-space and thesummation convention over repeated indices applies and will throughout.

The Green's tensor G_(αβ)(r,r′) satisfies the equation

∇×∇×G(r,r′)−k ₀ ² n ²(z)G(r,r′)=4πδ(r−r′)I  (23)

where I is the unit tensor. The Green's tensor must also satisfy theboundary conditions

{circumflex over (z)}×G(r,r′)|_(z=0) _(⁺) ={circumflex over(z)}×G(r,r′)|_(z=0) _(⁻) ,  (24)

{circumflex over (z)}×∇×G(r,r′)|_(z=0) _(⁺) ={circumflex over(z)}×∇×G(r,r′)|_(z=0) _(⁻) .  (25)

Making use of a planewave decomposition, it may be found that$\begin{matrix}{{{G_{\alpha\beta}( {r,r^{\prime}} )} = {\frac{i}{2\quad \pi}{\int\quad {\frac{^{2}q}{k_{z}(q)}{S_{\alpha \quad \gamma}^{- 1}(q)}{g_{\gamma \quad \delta}( {q,z} )}{S_{\delta\beta}(q)}{\exp \lbrack {{{ik}(q)} \cdot ( {r - r^{\prime \quad}} )} \rbrack}}}}},} & (26)\end{matrix}$

where it is assumed that z>z′>0. Here the matrix S(q) rotates k(q) intothe xz plane, or more explicitly $\begin{matrix}{{S(q)} = {{q}^{- 1}\begin{pmatrix}q_{x} & q_{y} & 0 \\{- q_{y}} & q_{x} & 0 \\0 & 0 & {q}\end{pmatrix}}} & (27) \\{{g_{xx} = {( \frac{k_{z}(q)}{k_{0}} )^{2}\{ {1 + {{R_{2}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \}}},} & (28) \\{{g_{yy} = {1 + {{R_{1}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}}},} & (29) \\{{g_{zz} = {( \frac{q}{k_{0}} )^{2}\{ {1 - {{R_{2}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \}}},} & (30) \\{{g_{zx} = {\frac{{- {q}}{k_{z}(q)}}{k_{0}^{2}}\{ {1 + {R_{2}(q){\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \}}},} & (31) \\{{g_{xz} = {\frac{{- {q}}{k_{z}(q)}}{k_{0}^{2}}\{ {1 - {{R_{2}(q)}{\exp \lbrack {2{{ik}_{z}(q)}z^{\prime}} \rbrack}}} \}}},} & (32)\end{matrix}$

all other elements of g being zero, R₁(q) is defined in equation (10),and R₂(q) is given by $\begin{matrix}{{R_{2}(q)} = {\frac{{k_{z}^{\prime}(q)} - {{nk}_{z}(q)}}{{k_{z}^{\prime}(q)} + {{nk}_{z}(q)}}.}} & (33)\end{matrix}$

We take the incident field to be an evanescent plane wave withpolarization E⁽⁰⁾

E _(α) ^(i)(r)=E _(α) ⁽⁰⁾exp [ik(q)·r],  (34)

where k(q) is the incident wavevector and k₀≦|q|≦nk₀. Within theaccuracy of the first Born approximation, we find that the scatteredfield is given by the expression

E _(α) ^(s)(r)=k ₀ ² ∫d ³ r′G _(αβ)(r,r′)E _(β)⁽⁰⁾exp[ik(q)·r′]η(r′).  (35)

In the far zone of the scatterer the Green's tensor assumes theasymptotic form $\begin{matrix}{{G_{\alpha\beta}( {r,r^{\prime}} )} \sim {{S_{\alpha \quad \gamma}^{\dagger}(q)}{g_{\gamma \quad \delta}( {q,z} )}{S_{{\delta\beta}\quad}(q)}\frac{\exp ( {{ik}_{0}r} )}{r}{\exp \lbrack {{- {{ik}(q)}} \cdot r^{\prime}} \rbrack}}} & (36)\end{matrix}$

where k(q) lies in the direction of r. Thus we find that the sccatteredfield becomes $\begin{matrix}{{{E_{\alpha}^{s}(r)} \sim {{A_{\alpha\beta}( {q_{1},q_{2}} )}E_{\beta}^{(0)}\frac{\exp ( {{ik}_{0}r} )}{r}}},} & (37)\end{matrix}$

where q₁ and q₂ are the transverse parts of the of the incident andoutgoing wavevectors, respectively. A_(αβ)(q₁,q₂) denotes the tensorscattering amplitude which is related to the susceptibility by

A _(αβ)(q ₁ ,q ₂)=∫d ³ rw _(αβ)(q ₂ ,z)exp{i[k(q ₁)−k(q₂)]·r}η(r),  (38)

where

w_(αβ)(q,z)=k ₀ ² S _(αγ) ^(†)(q)g _(γδ)(q,z)S _(δβ)(q).  (39)

Inversion of the above integral equation will be discussed in Section2.2.

2. Inverse Problem

The inverse problem consists of reconstructing the susceptibility frommeasurements of the scattering amplitude. To this end we will constructthe pseudoinverse solution of the integral equations (17) and (38). Webegin by briefly reviewing the singular value decomposition (SVD) oflinear operators on Hilbert spaces. Let A denote a linear operator withkernel A(x,y) which maps the Hilbert space H₁ into the Hilbert space H₂.By the SVD of A we mean a representation of the form $\begin{matrix}{{{A( {x,y} )} = {\sum\limits_{n}\quad {\sigma_{n}{g_{n}(x)}{f_{n}^{*}(y)}}}},} & (40)\end{matrix}$

where σ_(n) is the singular value associated with the singular functionsf_(n) and g_(n). The {f_(n)} and {g_(n)} are orthonormal bases of H₁ andH₂, respectively and are eigenfunctions with eigenvalues σ_(n) ² of thepositive self-adjoint operators A*A and AA*:

A*Af _(n)=σ_(n) ² f _(n)  (41)

AA*g _(n)=σ_(n) ² g _(n).  (42)

In addition, the f_(n) and g_(n) are related by

Af _(n)=σ_(n) g _(n)  (43)

A*g _(n)=σ_(n) f _(n).  (44)

The pseudoinverse solution to the equation Af=g is defined to be theminimizer of ∥Af−g∥ with smallest norm. This well-defined elementf⁺∈N(A)^(⊥) is unique and may be shown to be of the form f⁺=A⁺g, wherethe pseudoinverse operator A⁺ is given by A⁺=A*(AA*)⁻¹ and N(A)¹⁹⁵ isthe orthogonal complement of the null space of A. The SVD of A may beused to express A⁺ as $\begin{matrix}{{A^{+}( {x,y} )} = {\sum\limits_{n}\quad {\frac{1}{\sigma_{n}}{f_{n}(x)}{{g_{n}^{*}(y)}.}}}} & (45)\end{matrix}$

We now apply the foregoing to the inverse problem for scalar waves.

2.1 Scalar Case

The integral equation (17) may be rewritten in the form

A(q ₁ ,q ₂)=∫d ³ rK(q ₁ ,q ₂ ;r)η(r),  (46)

where the scattering operator K(q₁,q₂;r) is given by

K(q ₁ ,q ₂ ;r)=exp[i(q ₁ −q ₂)·ρ]κ(q ₁ ,q ₂ ;z),  (47)

with

κ(q ₁ ,q ₂ ;z)=k ₀ ²{1+R ₁(q ₂)exp[2ik _(z)(q ₂) z]}exp[i(k _(z)(q ₁)−k_(z)(q ₂))z]χ(q ₁ ,q ₂).  (48)

Here we have assumed that A(q₁,q₂) is measured for (q₁,q₂), in the dataset and have introduced a function χ(q₁,q₂) which is unity if (q₁,q₂) ∈and is zero otherwise. Note that specifies the available orientations ofthe incoming and outgoing wavevectors while χ(q₁,q₂) enforces theconstraints k₀≦|q₁|≦nk₀ and |q₂|≦k₀. In general, the wavevectors q₁ andq₂ vary continuously over all space. If, however, the transverse supportof η(r) is contained in the region [−L,L]×[−L,L], then q₁ and q₂ may betaken to be discrete and restricted to the set Λ={(n_(x)π/L,n_(yπ/L): n)_(x),n_(y)=0,±1, . . . }. Thus η(r) may be expressed as the Fourierseries $\begin{matrix}{{{\eta \quad ( {\rho,z} )} = {\sum\limits_{q \in \Lambda}{{c_{q}(z)}\quad {\exp ( {{iq} \cdot \rho} )}}}},} & (49)\end{matrix}$

where c_(q)(z) are appropriate coefficients.

In order to obtain the SVD of K(q₁,q₂;r) it will prove useful tointroduce the following identity: $\begin{matrix}{{{K( {q_{1},{q_{2};r}} )} = {\sum\limits_{Q \in \Lambda}{{\exp ( {{iQ} \cdot \rho} )}\quad \delta \quad ( {Q + q_{2} - q_{1}} ){\kappa ( {{Q + q_{2}},{q_{2};z}} )}}}},} & (50)\end{matrix}$

where δ denotes the kronecker delta. Using this result, we find that thematrix elements of the operator KK* are given by $\begin{matrix}{{{{KK}^{*}( {q_{1},{q_{2};q_{1}^{\prime}},q_{2}^{\prime}} )} = {\sum\limits_{Q \in \Lambda}{{M( {q_{2},{q_{2}^{\prime};Q}} )}\quad \delta \quad ( {Q + q_{2} - q_{1}} )\quad {\delta ( {Q + q_{2}^{\prime} - q_{1}^{\prime}} )}}}},} & (51)\end{matrix}$

where $\begin{matrix}{{{M( {q_{2},{q_{2}^{\prime};Q}} )} = {\int_{0}^{L}\quad {{z}\quad {\kappa ( {{Q + q_{2}},{q_{2};z}} )}{\kappa^{*}( {{Q + q_{2}^{\prime}},{q_{2}^{\prime};z}} )}}}},} & (52)\end{matrix}$

with L the range of η(r) in the {circumflex over (z)} direction. To findthe singular vectors g_(QQ′) of K which satisfy

KK*g _(QQ′)=σ_(QQ′) ² g _(QQ′);  (53)

we will make the ansatz

g _(QQ′)(q ₁ ,q ₂)=C _(Q′)(q ₂ ;Q)δ(Q+q ₂ −q ₁),  (54)

where Q,Q′∈Λ. Equation (51) now implies that $\begin{matrix}{{\sum\limits_{q^{\prime} \in \Lambda}\quad {{M( {q,{q^{\prime};Q}} )}{C_{Q^{\prime}}( {q^{\prime};Q} )}}} = {\sigma_{{QQ}^{\prime}}^{2}{{C_{Q^{\prime}}( {q;Q} )}.}}} & (55)\end{matrix}$

Thus C_(Q′)(q₂;Q) is an eigenvector of M(Q) labeled by Q′ witheigenvalue σ_(QQ′) ². Since M(Q) is self-adjoint, the C_(Q′)(q₂;Q) maybe taken to orthonormal. Next, the f_(QQ′) may be found fromK*g_(QQ′)=σ_(QQ′)f_(QQ′) and are given by $\begin{matrix}{{f_{{QQ}^{\prime}}(r)} = {\frac{1}{\sigma_{{QQ}^{\prime}}}{\sum\limits_{q \in \Lambda}{{\exp ( {{- {iQ}} \cdot \rho} )}\quad {\kappa^{*}( {{Q + q},{q;z}} )}{{C_{Q^{\prime}}^{*}( {q;Q} )}.}}}}} & (56)\end{matrix}$

It follows that the SVD of K(q₁,q₂;r) is given by the expression$\begin{matrix}{{K( {q_{1},{q_{2};r}} )} = {\sum\limits_{Q,Q^{\prime}}\quad {\sigma_{{QQ}^{\prime}}{f_{{QQ}^{\prime}}^{*}(r)}{{g_{{QQ}^{\prime}}( {q_{1},q_{2}} )}.}}}} & (57)\end{matrix}$

The SVD of equation (57) may now be used to obtain the pseudoinversesolution to the integral equation (46): $\begin{matrix}{{{\eta^{+}(r)} = {\sum\limits_{q_{1},\quad q_{2}}\quad {{K^{+}( {{r;q_{1}},q_{2}} )}{A( {q_{1},q_{2}} )}}}},} & (58)\end{matrix}$

where K⁺(r;q₁,q₂) is the pseudoinverse of K(q₁,q₂;r). Using the resultof equation (45), the pseudoinverse K⁺ may be seen to be given by$\begin{matrix}{{K^{+}( {{r;q_{1}},q_{2}} )} = {\sum\limits_{Q,Q^{\prime}}\quad {\frac{1}{\sigma_{{QQ}^{\prime}}}{f_{{QQ}^{\prime}}(r)}{{g_{{QQ}^{\prime}}^{*}( {q_{1},q_{2}} )}.}}}} & (59)\end{matrix}$

Substituting equations (54) and (56) into equation (59) and using thespectral decomposition $\begin{matrix}{{{\sum\limits_{Q^{\prime}}\quad {\frac{1}{\sigma_{{QQ}^{\prime}}^{2}}{C_{Q^{\prime}}( {q;Q} )}{C_{Q^{\prime}}^{*}( {q^{\prime};Q} )}}} = {M^{- 1}( {q,{q^{\prime};Q}} )}},} & (60)\end{matrix}$

where M⁻¹(q,q′;Q) is the qq′ matrix element of M⁻¹(Q) we obtain$\begin{matrix}{{{\eta^{+}(r)} = {\sum\limits_{q_{1},q_{2},q_{2}^{\prime}}\quad {\sum\limits_{Q}\quad {{\exp ( {{- }\quad {Q \cdot \rho}} )}{\delta ( {Q + q_{2} - q_{1}} )}{M^{- 1}( {q_{2},{q_{2}^{\prime};Q}} )} \times {\kappa^{*}( {{Q + q_{2}^{\prime}},{q_{2}^{\prime};}} )}{A( {q_{1},q_{2}} )}}}}},} & (61)\end{matrix}$

which is the inversion formula for scalar TIRT.

2.2 Vector Case

The integral equation (38) may be rewritten in the form

A _(αβ)(q ₁ ,q ₂)=∫d ³ rK _(αβ)(q ₁ ,q ₂ ;r)η(r),  (62)

where the scattering operator K_(αβ)(q₁,q₂;r) is given by

K _(αβ)(q ₁ ,q ₂ ;r)=exp[i(q ₁ −q ₂)·ρ]κ_(αβ)(q ₁ ,q ₂ ;z),  (63)

and

κ_(αβ)(q ₁ ,q ₂ ;z)=w _(αβ)(q ₂ ;z)exp[i(k _(z)(q ₁)−k _(z)(q ₂))z]χ(q ₁,q ₂).  (64)

As in the scalar case we have assumed that A(q₁, q₂) is measured for(q₁, q₂), in the data set and have introduced the appropriate blockingfunction χ(q₁, q₂). The vector integral equation (62) differs from itsscalar counterpart equation (46) only by a factor associated with thepolarization. Evidently, by measuring a fixed component of the scatteredfield for a particular incident direction we see that the scalarinversion formula (61) may be used to reconstruct η(r).

We now consider the SVD for the vector case. This is a simplegeneralization of the scalar case. Following the previous development wefind that the SVD of K_(αβ)(q₁,q₂;r) is of the form $\begin{matrix}{{K_{\alpha \quad \beta}( {q_{1},{q_{2};r}} )} = {\sum\limits_{Q,Q^{\prime}}\quad {\sigma_{{QQ}^{\prime}}{f_{{QQ}^{\prime}}^{*}(r)}{{g_{{QQ}^{\prime}}^{\alpha \quad \beta}( {q_{1},q_{2}} )}.}}}} & (65)\end{matrix}$

Here the singular functions are given by

g _(QQ′) ^(αβ)(q ₁ ,q ₂)=C _(Q′) ^(αβ)(q ₂ ;Q)δ(Q+q ₂ −q ₁),  (66)

$\begin{matrix}{{f_{{QQ}^{\prime}}(r)} = {\frac{1}{\sigma_{{QQ}^{\prime}}}{\sum\limits_{q \in \Lambda}\quad {{\exp ( {{- }\quad {Q \cdot \rho}} )}{\kappa^{*}( {{Q + q},{q;z}} )}{{C_{Q^{\prime}}^{*}( {q;Q} )}.}}}}} & (67)\end{matrix}$

The C_(Q′) ^(αβ)(q₂; Q) are eigenfuntions of M_(αβ) ^(α′β′)(q₂, q′₂; Q)with eigenvalues=σ_(QQ′) ² $\begin{matrix}{{{\sum\limits_{q^{\prime} \in \Lambda}\quad {{M_{\alpha \quad \beta}^{\alpha^{\prime}\beta^{\prime}}( {q,{q^{\prime};Q}} )}{C_{Q^{\prime}}^{\alpha^{\prime}\beta^{\prime}}( {q^{\prime};Q} )}}} = {\sigma_{{QQ}^{\prime}}^{2}{C_{Q^{\prime}}^{\alpha \quad \beta}( {q;Q} )}}},} & (68)\end{matrix}$

where $\begin{matrix}{\quad {{M_{\alpha \quad \beta}^{\alpha^{\prime}\beta^{\prime}}( {q_{2},{q_{2}^{\prime};Q}} )} = {\int_{0}^{L}\quad {{z}\quad {\kappa_{\alpha \quad \beta}( {{Q + q_{2}},{q_{2};z}} )}{{\kappa_{\alpha^{\prime}\beta^{\prime}}^{*}( {{Q + q_{2}^{\prime}},{q_{2}^{\prime};z}} )}.}}}}} & (69)\end{matrix}$

The pseudoinverse solution to the integral equation (62) is

given by $\begin{matrix}{{{\eta^{+}(r)} = {\sum\limits_{q_{1},q_{2}}\quad {{K_{\alpha \quad \beta}^{+}( {{r;q_{1}},q_{2}} )}{A_{\alpha \quad \beta}( {q_{1},q_{2}} )}}}},} & (70)\end{matrix}$

where $\begin{matrix}{{K_{\alpha \quad \beta}^{+}( {{r;q_{1}},q_{2}} )} = {\sum\limits_{Q,Q^{\prime}}\quad {\frac{1}{\sigma_{{QQ}^{\prime}}}{f_{{QQ}^{\prime}}(r)}{{{g_{{QQ}^{\prime}}^{\alpha\beta}}^{*}( {q_{1},q_{2}} )}.}}}} & (71)\end{matrix}$

More explicitly, we have $\begin{matrix}{{{\eta^{+}(r)} = {\sum\limits_{q_{1},q_{2},q_{2}^{\prime}}\quad {\sum\limits_{Q}\quad {{\exp ( {{- }\quad {Q \cdot \rho}} )}{\delta ( {Q + q_{2} - q_{1}} )}( {M^{- 1}(Q)} )_{\alpha \quad \beta}^{\alpha^{\prime}\beta^{\prime}}( {q_{2},q_{2}^{\prime}} ) \times {\kappa_{\alpha^{\prime}\beta^{\prime}}^{*}( {{Q + q_{2}^{\prime}},{q_{2}^{\prime};z}} )}{A_{\alpha \quad \beta}( {q_{1},q_{2}} )}}}}},} & (72)\end{matrix}$

which is the inversion formula for vector TIRT.

3. Regularization

In order to avoid numerical instabilty and set the resolution of thereconstructed image to be comensurate with the available data, the SVDinversion formulas must be regularized. In particular, we replace 1/a inthe inversion formulas (61) and (72) by R(σ) where R is a suitableregularizer. where R is a suitable regularizer. The effect ofregularization is to limit the contribution of the small singular valuesto the reconstruction. One way to do this is to simply cutoff all abelow some cutoff σ_(c). That is we set $\begin{matrix}{{{R\quad (\sigma)} = {\frac{1}{\sigma^{2}}\theta \quad ( {\sigma - \sigma_{c}} )}},} & (73)\end{matrix}$

θ denoting the usual Heavyside step function. Alternatively a smoothcutoff, such as that derived by the Tikhonov method, may be employed.

4. Example

To demonstrate the feasibility of the inversion using a scalarformalation, for simplicity, the reconstruction of η(r) has beenobtained for a collection of spherical scatterers. This collection isrepresentative of physical structures which may be imaged, such as anano-structured material; the collection presents the necessarydielectric contrast to effect direct reconstruction. The forward datawas calculated by considering the scattering of evanescent waves from ahomogeneous sphere including multiple scattering terms by means of apartial wave expansion. Consider a sphere of radius a centered at thepoint r₀ with refractive index n, n being related to the scatteringpotential or susceptibility by the expression n²=1+4πη. It may be foundthat $\begin{matrix}{{{A( {k_{1},k_{2}} )} = {{\exp \lbrack {{( {k_{1} - k_{2}} )} \cdot r_{0}} \rbrack}{\sum\limits_{l = 0}^{\infty}{( {{2l} + 1} )A_{l}{P_{l}( {{\hat{k}}_{1} \cdot {\hat{k}}_{2}} )}}}}},} & (74)\end{matrix}$

where A_(l) are the usual partial wave expansion coefficients, P_(l) arethe Legendre polynomials, and {circumflex over (k)}=k/{square root over(k·k*)}. Since evanescent waves are considered, the argument of theLegendre polynomials in equation (74) may exceed unity. The series maynonetheless be shown to be convergent due to the rapid decay of theA_(l) with increasing l.

The forward data was obtained for a collection of four spheres of radiusλ/20. All scatterers are present simultaneously in the forwardcomputation with inter-sphere scattering neglected. The incidentevanescent waves were taken to have transverse wave vectors q₁ on aCartesian grid with spacing k₀/4 such that k₀≦|q|≦nk₀, n being therefractive index of the prism used to generate the incident waves.Values of q₂ were also taken on the grid with spacing k₀/4 and |q₂|≦k₀,consistent with a measurement scheme in which the scattered field ismeasured in the far zone of the lower half-space. The spheres arearranged in two layers, one equatorial plane conincident with the z=λ/20plane, the other with the z=λ/4 plane. In each layer, one sphere wastaken to have susceptibility 4πη=0.44 (index n=1.2) and one sphere wastaken to have susceptibility 4πη=0.4+0.48i (index n=1.1+0.2i). In eachof the experiments, complex Gaussian noise of zero mean was added to thesignal at various levels. The experiments were performed for twodifferent prisms, namely, for FIG. 2A, the prism has an index of n=10,as might be encountered in the infrared range, and for FIG. 2B, an indexn=4, as might be encountered in the visible range. The data sets arethus composed of 4980×49 points for the higher-index prism, and 752×49data poins for the lower-index prism. The regularization parameter σ_(c)was taken to be 10⁻⁶ in computing the tomographic image at z λ/20 and10⁻⁵ at z=λ/4.

One can see from the reconstructions of FIGS. 2A and 2B that the realand imaginary parts of the susceptibility may be found separately andthat the reconstructions are subwavelength resolved. The resolutiondepends both on the size of the regularization parameter, whichindirectly sets the number of singular functions used in thereconstruction, and on the depth, a consequence of the fact that theprobe fields decay exponentially into the sample, resulting in loss ofhigh-frequency Fourier components of the spatial structure of thesusceptibility. The tomographs at the z=λ/20 layer are more highlyresolved for the higher-index prism than for the lower-index prism, butthere is little difference at the z=λ/4 layer.

5. System

As depicted in high-level block diagram form in FIG. 3, system 300 is atomography system for generating an image of an scatterer/object usingmeasurements of scattered waves emanating from an object in response towaves illuminating the object. In particular, object 100 is shown asbeing under investigation. System 300 is composed of: source 320 forprobing the object 100 through prism 105; data acquisition detector 330for detecting the scattering data corresponding to the scattered wavesfrom object 100 at one or more locations proximate to object 100;position controller 340 for controlling the locations of detectors 330and sources 320; and computer processor 350, having associated inputdevice 360 (e.g., a keyboard) and output device 370 (e.g., a graphicaldisplay terminal). Computer processor 350 has as its inputs positionalinformation from controller 340 and the measured scattering data fromdetector 330. Even though the scatterer is shown as being present inFIG. 3, actually two sets of measurements are obtained, namely, one setwith the scatterer removed, and another set with the scatterer present,to provide the necessary data for image reconstruction, as detailedabove.

Computer 350 stores a computer program which implements the directreconstruction algorithm; in particular, the stored program processesthe measured scattering data to produce the image of the object orobject under study using a prescribed mathematical algorithm. Thealgorithm is, generally, determined with reference to an integraloperator relating the scattering data to the forward scattering operatoras expressed by integral equation (17) or equation (38).

6. Flow Diagram

The methodology carried out by the present invention is set forth inhigh-level flow diagram 400 of FIG. 4 in terms of the illustrativesystem embodiment shown in FIG. 3. With reference to FIG. 4, theprocessing effected by block 410 enables source 320 and data acquisitiondetector 330 so as to measure the scattering data emanating fromscatterer 100 due to illuminating waves from source 320. Thesemeasurements are passed to computer processor 350 from data acquisitiondetector 330 via bus 331. Next, processing block 420 is invoked tocompute the kernel expressed by equation (17) or equation (38), whichmay for efficiency be pre-computed and stored. In turn, processing block430 is operated to execute the reconstruction algorithm set forth inequation (17) or equation (38), thereby determining the scatteringpotential η(r). Finally, as depicted by processing block 440, thereconstructed tomographic image corresponding to η(r) is provided tooutput device 370 in a form determined by the user; device 370 may be,for example, a display monitor or a more sophisticated three-dimensionaldisplay device.

Although the present invention have been shown and described in detailherein, those skilled in the art can readily devise many other variedembodiments that still incorporate these teachings. Thus, the previousdescription merely illustrates the principles of the invention. It willthus be appreciated that those with ordinary skill in the art will beable to devise various arrangements which, although not explicitlydescribed or shown herein, embody principles of the invention and areincluded within its spirit and scope. Furthermore, all examples andconditional language recited herein are principally intended expresslyto be only for pedagogical purposes to aid the reader in understandingthe principles of the invention and the concepts contributed by theinventor to furthering the art, and are to be construed as being withoutlimitation to such specifically recited examples and conditions.Moreover, all statements herein reciting principles, aspects, andembodiments of the invention, as well as specific examples thereof, areintended to encompass both structural and functional equivalentsthereof. Additionally, it is intended that such equivalents include bothcurrently know equivalents as well as equivalents developed in thefuture, that is, any elements developed that perform the function,regardless of structure.

In addition, it will be appreciated by those with ordinary skill in theart that the block diagrams herein represent conceptual views ofillustrative circuitry embodying the principles of the invention.

What is claimed is:
 1. A method for generating a tomographic image of anobject comprising: probing the object with incident evanescent waves,detecting waves scattered by the object, and reconstructing thetomographic image by executing a prescribed mathematical algorithm withreference to the incident waves and the scattered waves to generate thetomographic image with subwavelength resolution.
 2. The method asrecited in claim 1 wherein the detecting includes measuring scatteringdata from the object, said scattering data being related to a scatteringpotential of the object by an integral operator.
 3. The method asrecited in claim 2 wherein the reconstructing includes reconstructingthe tomographic image by executing the prescribed mathematicalalgorithm, determined with reference to the integral operator, on thescattering data, the prescribed mathematical algorithm further relatingthe scattering potential to the scattering data by another integraloperator.
 4. The method as recited in claim 1 wherein the probingincludes illuminating a prism to generate the incident waves.
 5. Amethod for generating a tomographic image of an object comprising:illuminating the object with incident evanescent waves, measuringscattering data from the object wherein the scattering data is relatedto the object by an integral operator, and reconstructing thetomographic image by executing a prescribed mathematical algorithm,determined with reference to the integral operator, on the scatteringdata to generate the tomographic image with subwavelength resolution. 6.The method as recited in claim 5 wherein the scattering data is relatedto a scattering potential of the object by the integral operator, andwherein the reconstructing includes reconstructing the tomographic imageby executing the prescribed mathematical algorithm, determined withreference to the integral operator, on the scattering data, theprescribed mathematical algorithm further relating the scatteringpotential to the scattering data by another integral operator.
 7. Themethod as recited in claim 5 wherein the illuminating includesilluminating a prism to generate the incident waves.
 8. A system forgenerating a tomographic image of an object comprising: a source forprobing the object with incident evanescent waves, a detector fordetecting waves scattered by the object, and a processor forreconstructing the tomographic image by executing a prescribedmathematical algorithm with reference to the incident waves and thescattered waves to generate the tomographic image with subwavelengthresolution.
 9. The system as recited in claim 8 wherein the detectorincludes means for measuring scattering data from the object, saidscattering data being related to a scattering potential of the object byan integral operator.
 10. The system as recited in claim 9 wherein theprocessor includes means for reconstructing the tomographic image byexecuting the prescribed mathematical algorithm, determined withreference to the integral operator, on the scattering data, theprescribed mathematical algorithm further relating the scatteringpotential to the scattering data by another integral operator.
 11. Asystem for generating a tomographic image of an object comprising: asource for illuminating the object with incident evanescent waves, adetector for measuring scattering data from the object wherein thescattering data is related to the object by an integral operator, and aprocessor for reconstructing the tomographic image by executing aprescribed mathematical algorithm, determined with reference to theintegral operator, on the scattering data to generate the tomographicimage with subwavelength resolution.
 12. The system as recited in claim11 wherein the scattering data is related to a scattering potential ofthe object by the integral operator, and wherein the processor includesmeans for reconstructing the tomographic image by executing theprescribed mathematical algorithm, determined with reference to theintegral operator, on the scattering data, the prescribed mathematicalalgorithm further relating the scattering potential to the scatteringdata by another integral operator.